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Abstract 

In this paper, we investigate moment methods from a general point of view using an 
operator notation. This theoretical approach lets us explore the moment closure problem 
in more detail. This gives rise to a new idea, proposed in |14l 115] . of how to improve the 
well-known Pjv approximations. We systematically develop a diffusive correction to the Pn 
equations from the operator formulation — the so-called Dn approximation. We validate the 
new approach with numerical examples in one and two dimensions. 

1 Introduction 

Developing simplified methods for the simulation of radiative transfer requires taking into 
account the physical situation that will be analyzed. There are two important limits: optically 
thick and optically thin media. In optically thin media there are very few particles that 
interact with the radiation. The distances that photons would typically travel before they 
are scattered or absorbed are therefore very long compared to the domain size. On the other 
hand, in optically thick regimes those distances are very short compared to the domain size. 




(a) Optically thin medium. 



(b) Optically thick medium. 



Figure 1 : Path of a photon in optically thin and optically thick media. 



The different regimes can be characterized by the mean free paths \ f° r scattering and 
e for absorption. Large mean free paths represent an optically thin regime, small mean free 
paths represent optically thick regimes. One problem in this characterization is, that many 
materials are optically thick in a specific frequency range and optically thin in other ranges. 
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Additionally, there is a transition regime between the two. This is the region of mean free 
paths between optically thin and thick media. In this regime, photons travel further than in 
optically thick situations, but not far enough as that the regime counts as optically thin. 

An example of an optically thin regime is radiation propagating in vacuum. Optically 
thick regimes can be found in glass cooling processes or combustion chambers. There are also 
situations where all these regimes play a role. For instance, during the reentry of a space craft 
into the atmosphere the regime goes from optically thin (space) through a transition (higher 
atmosphere) into optically thick (lower atmosphere). 

Due to large mean free paths in optically thin regimes it is possible to track the trace of 
many single photons until they leave the computational domain or undergo absorption. This 
is the basic idea of ray tracing and Monte Carlo methods. If there are very few scattering 
events, following the path of one photon is rather simple. These particle methods are often 
used in astrophysics where the physical conditions are usually in a way that these methods 
can be applied successfully. For examples see |18l 1231 1291 [30] . 

In optically thick regimes, following the path of a single photon is almost impossible 
because it undergoes too many scattering events before it reaches its destination (leaving 
the computational domain or being absorbed). Therefore, other methods are used. They 
are usually based on the diffusion approximation [20], e.g. the Simplified Pjv approximation 
[EHini ED HZ] or flux-limited diffusion [T2]. 

The transition regime is the region of optical depth that is located between optically thin 
and optically thick regimes. Methods that work well in optically thin media are computation- 
ally too expensive in these regimes. Methods that work well in optically thick media, on the 
other hand, give poor results for low order approximations and have high computational costs 
if one increases the order. Therefore, in these regimes new methods have to be developed. 

These new approximations have to recover traditional reduced models for small mean free 
paths. Further, in the transition regime, they have to be more accurate than the simplified 
models and should be solvable more efficiently than the full kinetic approaches. 

The field of use for transition regime models can be found in stand-alone solvers for 
problems that lie completely within the optically thick and transition regime. For problems 
where the order of the mean free path also covers the optically thin regime, the new approaches 
could be used in hybrid methods. 

In this work, our starting point will be moment methods. These are usually derived based 
on the assumption that in highly scattering materials the photon distribution is driven toward 
a local equilibrium Therefore, the radiative intensity distribution is almost isotropic at every 
point. If this is the case, instead of treating the full intensity distribution, we can restrict 
our analysis to quantities that are averages of the directional distribution function over all 
directions. These quantities are e.g. the spectral energy distribution, the spectral flux and 
the spectral pressure. Averages of products of intensity distribution with directional test 
functions are called moments of the intensity function. Often, one is only interested in these 
averaged quantities. 

There is a variety of moment closures. One of the first was the Pn closure, which was 
developed by Chandrasekhar Other approaches include the minimum entropy closure 
E E 1171 1121 113] . A recent approach applies methods from the study of dynamical systems to 
moment closure [9l 122]. 

In Section [2] we introduce the main concept of moment methods in an operator notation. 
This theoretical approach lets us explore the moment closure problem in more detail. This 
investigation gives rise to a new approach, proposed in |14l 115] . of how to improve the well- 
known Pjv approximations (Section [3}. Using the operator formulation, we systematically 
develop a diffusive correction to the Pjv equations — the so-called Dn approximation. We do 
not treat boundaries here; they are considered in [8]. The partial differential equations behind 
the operator notation are developed for a simple example in Section [3] before we develop the 
general Pjv and Dn equations in Sections[S]and|S]respectively. Numerical examples in one and 
two dimensions are then shown in Section [7] We find that the solution of the Dn equations 
is at least as accurate as the solution of the Pjv equation of two orders higher. 
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2 Moment Models and Deviation Decomposition 



We consider radiation in a spatial domain X with boundary dX whose intensity I(t, x, Q) 
at time t > 0, position x £ X, and direction Q, € S 2 is governed by the frequency-averaged 
radiative transfer equation (RTE) 

-d t I(t, x, fi) + ft ■ V* 7(t, a;, ft) + (o-(x) + K{x))I(t, x, fi) 



ct(x) 



/ Q.n')/(£, a;, fiV n ' + «(a:)-B(T(a;)) + Q(t, a;,Sl). (2.1) 

Js 2 



Here cr(x) and k(x) are the scattering and absorption coefficients, &(x, S7 • S7' ) > is the scat- 
tering redistribution function, B(T(x)) > is the blackbody emission intensity at temperature 
T(x) > 0, and Q(t, x, Q) is the emission due to other sources. The scattering redistribution 
function satisfies the normalization 

$(x,Sl • S7')dS7 = 1. (2.2) 

The fact that a(x), k(x), and B(T(x)) are independent of SI, while &(x,Q ■ S7') depends on 
Q ■ S7' is consistent with a stationary, isotropic background medium. The fact that these 
functions are independent of t means that the heat capacity of this medium is large. 
We will express the different parts of equation (|2.1|l in terms of operators. 

Definition 2.1. We define the operators 

(Ai)(t,x,n) = n-v x i(t,x,n), (2.3a) 

(SI)(t,x,fl) = f $(x,n- n')I(t,x,Q')dQ.' , (2.3b) 

4tt J S 2 

(!CI)(t,x,n) = (n(x) + a{x))I(t,x,Q) - (SI)(t,x,Q) , (2.3c) 

(CI)(t,x,Q.) = (A+fC)I(t,x,n) , (2.3d) 

Q(t, x, fi) = k(x)B{T) + Q(t, x, SI) . (2.3e) 

Here A is advection, S is scattering, K, is total interaction due to scattering and absorption, 
and Q(t,x,Q.) is total emission. Using these operators the RTE reads 

-d t I(t, x, S7) + CI(t, x, Q) = Q{t,x,Q). (2.4) 
c 

We impose homogeneous boundary conditions. Let 

V = dX x S 2 , and = {(ar, fi) € T : ±n(x) • SI > 0} , (2.5) 
where n is the outward unit normal vector. Appropriate boundary conditions are 

l(t,x,n)=0 Vf > and (x,fi) £ T~ . (2.6) 
Together with the initial condition 

I(p,x,Q) =I (x,n), (2.7) 

we have a well-posed problem. Under certain (physically reasonable) assumptions on the 
scattering and absorption coefficients a and re, and for Q £ L 2 (]0, t\ [xX x S 2 , K), Q(i,a;,n) > 
there exists a unique solution / € {/ £ D t : 7 = on T~} (cf. 0), where 

D,cL 2 (]0,ti[xXxS 2 ,K). (2.8) 

Moment methods have a long history [4] [6]. Nevertheless they are still used for solving 
radiative or neutron transfer problems in situations where computational time is of concern. 
The main idea of moment methods is to derive an approximation to the radiative intensity 
distribution with respect to its directional moments. This relation is used to find an expression 
for the closure relation. There are several ways to find such approximations. The Pjv approach 
expresses the radiative intensity distribution as a series expansion of spherical harmonics. 
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The minimum entropy approaches use an expression with an exponential function. But, 
independent of how these methods approximate the intensity, all of them have in common that 
the unknown coefficients in their expansion are somehow related to the directional moments 
of the intensity function. 

Moments are directional averages of the intensity distribution multiplied with a test func- 
tion that depends on the direction SI. As test functions one can choose between several 
possibilities. We will use spherical harmonics, denoted by Y k , cf. Appendix [A] There are 
several reasons for choosing these functions. First, they form an orthonormal basis of the 
function space H m (5' 2 ,C) and therefore, after transformation also of L 2 (S' 2 ,R). Hence, they 
can be used for a complete description of the dependence of the intensity distribution on 
the direction SI. In other words, each intensity distribution can be represented by the series 
expansion 

oo I 

I(t,x,Sl) = J2Y, , (2.9) 

!=0 k=-l 

with moments 



Ii(t,x)= [ Y k (Sl)I(t,x,Sl)dSl . (2.10) 
Js 2 



To deal with spherical harmonics and the related moments we define 
Definition 2.2. Moments of the radiative intensity are generated by the operator 

mf :D t -> L 2 (]0,ti[xX,C) 



'■ 



I(t,x,Q,)>-> Y k (S2)I(t,x,Sl)dSl , 



(2.11) 



and we write for the moment of order I and degree k 

I?(t,x)=m? (I(t,x,Sl)) . (2.12) 

We allow the moments to be complex valued. A real-valued approximation to the radiative 
intensity is obtained by taking the real part of these equations. The moments depend on time 
and space. For convenience we neglect this in the notation if it is clear to what we are referring 
and write if = l k (t,x). 

A vector I of all moments belongs to 

= {i = (..., If ,.. .) T : I e No, k G {-l, ...,/}} C I 2 (L 2 (]0, ii[xJf, C)) , (2.13) 



where I 2 denotes the space of all square summable sequences. We introduce 
Definition 2.3. We define the "Intensity to Moment" operator 

M : D t M t 

(2.14) 

I(t,x,Sl)^l(t,x) . 

The inverse transformation is given by the "Moment to Intensity" or "Expansion" operator 

£ : Mt -> D f 

^ ^ k » (2-15) 

1=0 fc=-i 

By their construction the operators M and £ are linear, bounded and continuous. Fur- 
thermore, it is easy to see that both operators are bijective. 

So far we have replaced the unknown dependence in SI by infinitely many unknown mo- 
ments. This does not help us to solve the RTE. A usual approach to overcome this problem 
is to assume that finitely many moments are sufficient to describe the intensity function. 
This reduces the amount of unknowns to a finite number and the problem can be handled 
much more easily. Assuming that only the moments up to order iV are relevant gives an 
approximation Ijsr(t,x, Si) to I(t,x,Sl) 

N I 

I(t,x,Sl)&I N (t,x,Sl) = J2 J2 I?(*,aO*i*(fi)- (2-16) 
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and it holds 



lim I N (t,x,Q) = I(t,x,Q.). (2.17) 

N— too 



The finite set of moments can be represented by the vector 

Ijv = (Io,I ■ ■ • , Ijv \ljv) T (2-18) 
and we define the set of restricted vectors of moments as 

Mf = jliv G (L 2 (]0,t 1 [xX,C)) (N+lf } • (2.19) 

Note that Mf is isomorphic to a subspace of Mi. 

We restrict ourself to approximations of the radiative intensity of odd orders. There 
are several reasons for this. First of all, even order approximations do not contain more 
information than odd order approaches. Therefore, they only introduce more moments and are 
computationally more expensive without giving any advantage. A second point for choosing 
just odd order approaches is given in [6| Chapter 10, § 3.2]. There it is shown, that boundary 
conditions for even order approximations are much less accurate than for odd order models. 

Analogous to Definition 12.31 we define 

Definition 2.4. The "restricted Intensity to Moment" operator is 

Mn :D t -cMf , , 

(2.20) 

l(t,x,0.) i— > Ijv(t,a;) . 

The inverse transformation is given by the "restricted Moment to Intensity" operator 

E N : Mf -> Of C D t 

N 1 (2 21) 

I N (t,x)^Yl X Ii(t,x)Yi k (ty = lN(t,x,n) . 



with 



= <In G D t : I N (t,x,Q) =$3 X) lUt,x)Yi k (n)\ . (2.22) 

I 1=0 k=-l J 

The only difference between the operators £ and Sn is the restriction on the domain and 
the range. The restriction of Range(£ jv) on Df ensures that the injectivity is inherited from 
£. Therefore, £n is still bijective. Df is the subspace of Dt that contains only those intensity 
functions, which can be represented by moments up to order N. Due to the bijectivity of 
£n, working with either the set of moments up to order iV or the approximated intensity 
distribution lN(t,x,Q) is equivalent. 

Lemma 2.5. The combined operator 

Vn ■ Dt -> D t 

I(t,x,Q) ^ £ N M N I(t,x,n) ( ' ' 

is a projection. 

Proof. We have to show that Pj V 7(n) = T-'jv/(fi) = ijv(fi) holds. Writing down the intensity 
function as series expansion, applying the operators as defined in Definition 12. 41 and addition- 
ally using the orthonormality property of the spherical harmonics gives the result. □ 

By using Vn, we can define the projection onto the orthogonal complement Dj _L Dj 
(with D t = Df Df) by Vn = Id -Vn- This gives rise to 

Definition 2.6. The radiative intensity can be decomposed into a component that can be 
described by finitely many moments 7jv(f2) and a deviation 

7jv(fi) =VnI{V) . (2.24) 

This splitting is called deviation decomposition. 
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We call In (fi) the deviation because it is the difference between the full radiative intensity 
7(fl) and the component Jjv(fi) that can be represented by finitely many moments: 

j(n) = r N i(n) + (id--Piv)i(n) 

= Jjsr(fi)+Pi\r/(fi) (2.25) 

Lemma 2.7. The RTE can be decomposed into an equivalent coupled system of finitely many 
moment equations and one deviation equation 

-d t I N + MnCSnIn + MnCIn = Qn , (2.26a) 

c 

-dJ N (n) + PnjCSnIn + VnCIn = Qn ■ (2.26b) 

c 

Proof. We start by decomposing the intensity distribution and source terms into one compo- 
nent that depends on a finite set of moments and second part that is the deviation 

/ = In + In = SnIn + In ■ (2.27) 
Q = Qn + Qn = SnQn + Qn . (2.28) 

By Qjv we denote the set of moments generated from the source term and Qn represents the 
deviation part 

Qn = {m° Q,m7 1 Q,...,m%Q,) T and Qn=TnQ. (2.29) 
Due to the linearity of the operators, this leads to 

-dt (SnIn + In) + CEnIn + Un = SnQn + Qn . (2.30) 

Applying the operator A4 jv to this equation gives (|2,26a[l . We have the orthogonality relations 
Ijv_L/jv and Qjv-LQjv that lead to 

M n (In)=0, Mn(Qn) = 0- (2.31) 

Using the operator Vn with (|2.30[) leads to ()2.26b|) . This is true since 

Vn{£n1n) = 0, Pn(SnQn) = . (2.32) 

□ 

We transformed the problem of solving the RTE from solving one equation in six dimen- 
sions to a different problem with a system of equations for the moments in four dimensions 
and additionally one coupled deviation equation that still is six dimensional. So far we have 
not gained anything. But if we could find a good approximation to the deviation, system 
(|2.26[) would simplify to just the moment equations where the dependence on the deviation 
could be treated explicitly. 



3 Closure Approximations 

The classical Pjv approximation is obtained by setting the deviation to zero 

In = . (3.1) 

This is equivalent to assuming that the radiative intensity distribution can be written as a 
finite sum of spherical harmonics. It is the simplest closure approximation one can make. Of 
course, in reality this is usually not true and thus, this assumption defines the limits of the 
Pjv approach. The Pn equations in operator notation are 

-d t I N + MnC£ n I n = Q N - (3.2) 
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In this work we derive a better approximation of In from the deviation equation. Starting 
from (|2.26b|) . we first assume that we can drop the time derivative of the deviation, whereby 

P N C£ N l N + P N CI N = Q N ■ (3.3) 

Then using the definition of the operator C = A + IC and assuming the invariance of IC under 
the projection Vn we obtain 

(VnA + tdjI N = Q N - TnCSnIn ■ (3.4) 

The invariance assumption is justified if the scattering kernel can be expanded in terms of 
spherical harmonics. For the second component on the right-hand side then holds 

VnjCEnT-n = VnAEnT-n + 1CPn£nIn 

= VnA£nIn +1C{M~Vn)£nIn (3.5) 
= VnASnIn + IC (In — In) = VnAS jvIjv ■ 

We thereby obtain 

In = {PnA + /C) (Qn - VnAEnIn) (3.6) 
which is a formal expression for the deviation. Of course, computing the inverse operator 
(PnA + }C\ is still not easier than solving the original transport equation. 
Starting from (|3.6jl and using a reformulation gives 

Mfi) = (id- (-IC-'-pNA^y 1 IC' 1 (Qn-VnA£nIn) . (3.7) 

Recall that we are interested in methods for the transition regime. Then the collisional physics 
are more important than the free transport of the photons. Hence, we can assume that the 
VnA component in (|3.6[) which represents free transport is significantly smaller than the IC 
component which describes absorption and scattering. We therefore formally use Neumann's 
series to obtain 

oo 

In (SI) = ^{-K.- 1 f N Ay IC' 1 (Qw(fi) -TnASnIn) ■ (3.8) 

3=0 

Of course, the operator A is not bounded and in general this series will not converge. Nev- 
ertheless, truncating the expansion after terms of some order gives an approximation to the 
deviation that can be used in the system of moment equations (|2.26a|l to improve the re- 
sults compared to the Pjv method. In this work we will deal with the approximation that is 
obtained by taking only the first term of (|3.8|l . This leads to 



In = IC' 1 (Qn - PnASnIn) ■ (3.9) 



Additionally, we remark that due to the orthogonality of the two projections VnI and VnI 
we have 

MnK.{kt x (Qn - VnASnIn)) =0- (3.10) 

Using the deviation approximation (|3.9p in (|2.26a|l finally leads to the Dn equations in oper- 
ator notation 

-d t I N + MnCSnIn + MnA^IC' 1 (Qn - TnASnIn)) =Qn- (3.11) 

Instead of computing the inverse of the operator (t>nA + IC\ we now only have to express 

the inverse of the combined absorption and scattering operator IC -1 . But this can be done in 
a straightforward way, cf. [B] 

The approximation of the deviation can be extended to higher orders. Truncating the 
series after terms of order zero gives an additional term with second derivatives in the moment 
equations. This is obvious since the operator A is applied twice. Using truncations after terms 
of higher order leads to deviations of higher orders and therefore, makes the equations much 
more complicated. 
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4 An Example 



In the next section, we are going to develop explicit expressions for the introduced operators 
for radiative transfer in 3D. But before we do this, to get an understanding how all these 
operators act on the equations, we take a closer look on a rather simple problem. We assume 
a one- dimensional slab geometry, i.e. the analyzed radiation field is homogeneous in two 
directions x\ and xi and also rotationally invariant with respect to the axis of propagation. 
Then the angular dependence can be expressed in one variable (i £ [—1,1]. For moments of 
the radiative intensity it holds 

Js 2 Jo J-i 



4-7T (n + m] 



a 



-i 



But the first integral is zero for every m^O and the set of relevant moments simplifies to 

for m = , 



(4.2) 

otherwise . 

Due to the homogeneous setup, all derivatives in direction of xi and X2 vanish and the radiative 
transfer equation becomes 

-d t I(t,x,fj)+txd X3 I(t,x,tJ,) + (a+K)I(t,x,fx)-(SI)(t,x,tx) = Q(t,x,/j) x <E X C K . (4.3) 

Moment equations can be generated by applying the operator to this equation: 

-ftl? + d X3 (fcs(0, 0I?+i + W)I?-i) + = Q? , (4.4) 
c 

where 

km ^7w=i' h3{0 ' l) = v(«+t)(a+3) ' and dl= a+K ~r i > (4 - 5) 

with o i being the moments of the scattering kernel, cf. Appendix [B] If in addition we assume 
an isotropic source, all moments of the source term of order unequal to zero vanish (Q; = 
for I > 0). It can be easily checked, that for the Pn approach this leads to the following set 
of equations 

-Otto + MO, 0)^ 3 I? + a T° = Q[J , (4.6a) 
c 

i<9 t I? + d X3 ha(0, 0^ 3 I?+i + h(0, l)d X3 tU + ffil? = , I e {2, . . . , N - 1} (4.6b) 

iftlS, + I 8 (0, JV)e- 3 Ijv-i + ^1% = • (4.6c) 

c 

For the Djv approach, we have to evaluate the expression 

M N A (k,- 1 (q n - V N M N I N y) . (4.7) 
Assuming an isotropic source Q gives Qjv(n) = 0. The moment to intensity operator becomes 

N 

s N i N = J2 y i 0i i > ( 4 - 8 ) 

1=0 

and by using the projection property of Vn we get 

N 

VnAEnIn = VnJ2 fW(n)0* 3 i° 

N (4 9) 

= Vn (fts(0,I)5l+i(n) + 's(0, 0*i-i(O)) 1? 

= / l3 (o,iV)y J ? +1 (n)9 !C 3iS, . 
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m 



(4.10) 



Applying A and K, 1 to this expression yields 

vCAK^VnAS nIn = m° n (d X3 (^—h 3 (0,N)n 3 Y N+1 (tyd X3 I N 

° (d X3 (^—h 3 (0, N)(h 3 (0, N + l)Y& +a (n) 

+ / 8 (o,Jv + i)y^(n))a,,i§ r 
(^7^(0, iv); 3 (o,iv + i)a.3i^) for n = jv, 

otherwise . 

From thes calculations we see that the Dn equations differ from the Pjv equations (|4.6[) only 
in the the equations for the moment of order N. The Dn system finally reads 

-d t I° + h 3 {0, Oj&jlS +aoI° = Qo , (4.11a) 

c 

i&i? + d X3 h 3 (o,i)d X3 i} +1 + i 3 (o,i)d X3 iti + = o , (4.11b) 

iftlS, + i s (0, iV)^ 3 I^-i - ^ 3 f-J— ^(0, JV)/ 3 (0, N + l)d X3 I° N ) + a N I° N = . (4.11c) 

C \<JN+1 J 

The correction term of the Dn equation is of diffusive nature and thus adds a stabilizing 
component to the Pjv equations. 

Remark 4.1. The additional term also can be interpreted in a different way. If we take the 
moment equation of order N + 1, neglect moments of order N + 2 and the time derivative 
and solve this equation for the moment we get 

I N+1 =-^—l 3 (0,N + l)d X3 I N . (4.12) 

Inserting this term as approximation for 1 N+1 into the equation for the moment 1^ gives 
exactly the equation for the moment of order N in the Dn equations. Therefore, at least in 
ID there is a simple way how the new model equations can be derived. 



5 Explicit Operators for Pn 

In the previous sections, we developed an operator approach to solve the RTE by moment 
methods. Now we take a more detailed look at these operators and analyze their structure. 
Also the results presented here are relevant to develop numerical methods for solving the RTE 
with the help of Pn and Dn equations. Most of the notation used here is introduced in the 
Appendix. 

In this section we will assume that the source term Q is isotropic and thus does not depend 
on the direction of the radiation fi. Then, due to the orthogonality of the spherical harmonics, 
all directional moments of Q(t,x) of order equal or higher than one vanish and we get 

oo ( 

Q(t,x) YfWmtQ^x) = -±=q%(t,x) . (5.1) 

tlkt^i V4^ 

For the vector of moments of the source and its deviation this leads to 

Q N (t,x) = (y^(K,B(T) + Q(t,x)),0,0,...y and Q N (t,x,Q)=0. (5.2) 
Next, we will analyze the expression 

MnCSnIn = MnASnIn + MnK,£nIn (5.3) 

which appears in the Pn and Djv approaches. The analysis will be performed separately 
for the transport (A) and the scattering/absorption (IC) component. We will start with the 
transport term A4nA£ jvIjv. 
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By analyzing one single moment equation of order n and degree m we get 

\r=l / \r=l \!=0 fc=-i 

3 / N I 

= EM E E (nca-y« fc (n)) if 



AT i 

EI 

r=l \i=0 fc=-! 

The inner sum can be written as a scalar product 

jv i 



(5.4) 



Y, E (nCO.r ; fc (fi)) If = (m^^Yjv, Ijv) = (nC^Yjvf Ijv , (5.5) 



i=0 fe=-Z 



where Yjv denotes the vector of spherical harmonics as introduced in Definition I A. 2 1 

Expressing the first component of this scalar product with the help of relation (|A,19[) and 
orthogonality relations for spherical harmonics gives for one component of the vector 



= 7- ((e|I, 3) ) T i:m™Y w + (e^ 1 ) r ^ +1 lCY Jf+1 j (5.6) 

( ( N \ T t N N ( JV+l\ T r >JV+l AT+1 1 

In our situation holds j, n £ {0, . . . , N} and therefore the unit vectors e^+j 1 and e^ 1 ^ are 

always zero in the last 2N + 3 components. Hence we can reduce the dimension of the last 
term without losing any information and end up with 

nC^Yjv = 7r (Ll + #£) e^, n) . (5.7) 

Finally we get for one moment equation 

3 



r=l 
3 



= E( (>(^ + C)e^,„)J (5.8) 
= E 7r ((-^0 +(y*r) ) d Xr I N ■ 

r=l ^ ' 

Taking advantage of the symmetry properties provided in Lemma lA.31 gives rise to the defi- 
nitions 

C X1 = L X1 + U X1 , C X2 = — L X2 — U X2 , C X3 = L X3 + U X3 . (5-9) 
For the set of all moment equations holds 



3 . 

MnASnIn = r )rC Xr d Xr lN = —C x1 8 x1 In + —C X2 X2 1n + C X3 8 X3 In ■ (5.10) 



Let us now analyze the scattering and absorption component A4 nICEnIn- As defined in (|2.3p 
it decomposes into 

ICS jvIat = (ft + er) £ atIjv — SSnIn (5-H) 
and from Appendix [B] (especially from ()B.7|) ) we get for one single moment and I < N 

W$K£nIn = ailf. (5.12) 

For the set of all moment equations we end up with 

MnICSnIn = Eivljv (5.13) 
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where we define 

/So 



V 



(5.14) 



with diagonal submatrices 

R 2i+i x R 2i+i B = ~, w ^ ^ 

Finally, the full set of moment equations for a Pjv approximation reads 
11 i 

—dtlN + TrC xl d xl lN + -^CxodxilN + C X3 d X3 lN + EjvIjv = Qjv (5.16) 
c Z Z 

with the vector of moments of the source term Qjv. This representation is equivalent to 
the operator notation given in (|3.2p but the operators are made explicit by their matrix 
representation. 

Proposition 5.1. The system, matrices C X1 and C X3 are symmetric. C X2 is skew symmetric. 
This means that the Pn equations are hyperbolic. 



Proof. Due to the symmetry results from Lemma lA.3l and the structure of the system matrices 
given in (|5.9p the result is obvious. □ 

6 Explicit Diffusion Operators for Dn 

The Dn equations (|3.11[1 have one additional component that was not already treated for the 
Pjv equations in Section [5] 

M N A (AT 1 (Qjv - VnMnIn)) ■ (6.1) 

In this section we investigate the properties of this term and develop a representation that 
fits into the matrix framework developed for the Pjv equations in (|5.16[) . 

Dealing only with isotropic source terms results in vanishing deviations Qn of the source 
term (see (|5.2[l ) and simplifies the problem to analyzing 

-MnAJC^VnMnIn ■ (6.2) 

We start with 



P n M n In=VnAJ2 Y, Yi k (tyli =-PnA(Y n ,I n ) 

(6-3) 

= VN^2(n s Y N ,d Xs i N ) . 



Formally the projection Vn is only defined as an operator that acts on functions from Dt into 
Dt. When working with vectors of functions we just apply the projection to every component. 

Using the decomposition of the orthogonal projection into Vn = Id— Vn and taking a 
component wise look at Vn^sY n leads us to 

V N n s Y; = V N (y. ( (e^Y LlY N + (e?+}) T U» +1 Y n+ i 



= 7, ((«#,,■)) T L^VnYn + (e^Y U?+ 1 V N Y N+ " 

Applying the projection operator Vn on spherical harmonics up to order N is an identity 
operation and it holds VnY™ = Y™ for n < N. If we apply Vn on spherical harmonic of 
order larger than N the result is VnY™ = for n > N. 

VnYn+i = (Y °, Yr\..., Y^-\Y^, 0, . . . , ) T = Y^+i ■ (6.5) 

2JV+3 zeros 
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Using these results for the orthogonal projection leads to 
Vn^sY; = (Id-Pjv^Yj 

" (Yjv- Yjvi 



7s (Yjv — Yjv) + (e^f U^ +1 (Yjv + x - Y^jr)) (6 



6) 



-y ( e N+1 \ T fj N+1 (() ny-"- 1 y n+i\ 



By using the fact that j < N and the special structure of the matrices U^ s +1 (see (|A.16[) ) we 
finally obtain for the full vector of spherical harmonics Yjv 



■pjv^sYjv = r y s Z^ lt Ui v s +1 Y N+ i 



with 



/ 











Id 
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(6.7) 



(6.8) 



\ o ■■■ o 

It holds Z? ut G R (JV+1)2 x R (JV+2)2 and the block with the identity matrix is located in the 
rows (iV 2 + 1) till (N + l) 2 and the columns (N 2 + 1) till (N + l) 2 . Finally we get 



J / a. \T 

VnAEnIn — ^2ls (z^ ut U^ s +1 Y jv+ij d Xs lN 

s=l 

= (Yjv +1 ) T E7, (t>^ +1 ) T (z^) T d Xs I> 



(6.9) 



Here we see, that only the moments of order N influence the additional term in the Djv 
equations. 

Applying the inverse of the combined scattering and absorption operator /C _1 and using 
results from Appendix [Bl gives 



By considering the complete term from (|6.2|) . for one single moment equation holds 

= IX (777 (^v + T ) I> (^ +1 ) T (z^) T d Xa h 

r=l \ s—1 



(6.10) 



(6.11) 



where we have to find an expression for m™Q r (Yjv+i) . Again, we start be analyzing on 
single component of the vector with n £ {0, . . . , N}, m £ {— n, . . . , n}, j (E {0, . . . , N + 1} and 
i G {-3, ■ ■ ■ J} 



JV+1\ fN+1 



L x „ Yjv+i + ( e 



JV+2 \ f T N+2 
=(i,3) 



N+2 



((e^fC 1 (hCYjv +1 ) + (e-+ 2 ) T U^ 2 (nCY„ +2 )) (6.12) 



Tr e 



N+1\ T fJV+1 JV+1 / JV+2\ T frJV+2 JV+2 



'(m,n) ^ V ^J) 



(m.n) 



Due to the orthogonality of the spherical harmonics the expressions m™Yjv+i and m^Yjv+2 
lead to unit vectors e^ 1 "^ and e^7^)- But since n < N the component that is one is always 
located in the first (N + l) 2 components. For the same reason (j < N + 1) the last 2jY + 4 
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components of e^+ 2 are always zero and we therefore can neglect the last 2N + 4 rows and 
columns of U£l +2 . Due to the structure of these matrices this is equivalent to writing 

m n -.A j AT+l\ T / f N+l , r~ T N + l\ JV+1 ( R -I o\ 

in„ n r Y 3 = 7r (e^J) ) [L X J + U Xt 7 ) e (n £ n) . (6.13) 
For a complete vector of spherical harmonics we obtain 

(6.14) 
(6.15) 

(6.16) 

(6.17) 
(6.18) 

(6.19) 



„™n v I fN+1 , t "tN+1\ N+l 

m„ n r Y N+1 = j r [L x ; + U Xr ^ ) e ( J n) . 
For the transposed expression that is relevant in (|6.11|l we end up with 

„"in l^r \T ( N+l \ T / f N+l , r"rN+l\ T 

tn n fi r (Yjv+i) = 7r (e (m *; n) J \L X ? +U mt T ) ■ 
If we use the fact L h Xr L% t = and U Xr U Xs = for k £ N we get 
m™AK. 1 VnA£nIn 

r — l \ s — 1 



or if we use 



with 



(l<* +1 ) 2 ) 2 3D r 



l x K => -^restrict 



/ 1 

V 



we end up with 



MnAK^VnMnIn = ~y^dx r ( rJ— y^V r , s d Xe lN | 



for the set of all moment equations. 

Finally, the moment equations for the Dn approach read 
11 i 

—dtT-N + — C x1 8 x1 In + -C X2 d X2 lN + Cx 3 d X3 lN 



ON+1 



(6.20) 



Due to the special structure of the matrices Ux a +1 and from (|A.16[) we get for the 

product 

, T 



(Li +i )\ux N s +i ) T =(K +i Li +i y 



For the complete matrix T> r ,s we obtain 



\ 



(6.21) 



({Ul r Ll s f 



f 



(6.22) 



V 



7r7s (U» +1 L» +1 ) T J 



From this matrix we see, that the Dn equations differ from the usual Pjv equations only in 
the equations for the moments of order N. 
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7 Numerical Results 



In this section we present some numerical simulations of two radiative transport problems. 
We compare the Dn method with the standard Pjv method. We first investigate a one- 
dimensional test problem that can be solved analytically. The second test problem compares 
the two methods in a strongly inhomogeneous two-dimensional medium. 



7.1 Pjv an d -Djv models vs. benchmark solution 

In this numerical test we compare the Pn and Dn approximations to an analytic benchmark 
solution. The details regarding the benchmark solution can be found in |26| . 

The physical setting is initially a cold, homogeneous, and infinite isotropically scattering 
medium. An internal slab radiation source is switched on at time t — and off at t — t cn d. 
Because this problem has slab symmetry, it can be described by the ID radiative transfer and 
energy transfer equations 

-d t I + (J0 X I + (a + k)I -SI= KaT 4 + Q , (7.1a) 

QCv iw = K (/ ~ 2aTi ) ■ (7 - lb) 

Assuming that the coefficients of absorption and scattering are constant and that 

c v =T 3 , (7.2) 

the usually nonlinear system becomes linear in the variables I and T 4 . For convenience we 
set c = 1 and g = 1. 

We impose the boundary conditions 

lim I(t,x,iMn) = , lim T(t,x) = 0, (7.3) 

and the initial conditions 

I(t = 0, x, n) = , T(t = 0,x) = 0. (7.4) 

A uniform isotropic radiation source is turned on in the slab [— xq,xq] over the time interval 
t £ [0, tend]- It can be described by 

Q(M,„) = M f ° r * e [-*„,*<,] and te[0,t enA ], 
I otherwise . 

Several benchmark solutions were provided in [26] for various values of the absorption and 
scattering coefficients. For our comparisons we will use the one with 

K = 1 , (7 = 0, XQ = 0.5 , tend = 10 . (7-6) 

The comparison of this benchmark solution with the Pjv and Djv solutions is given in Figure[5] 
The results have been obtained with a kinetic scheme for the transport part of the equation 
and a standard finite differences discretization of the difusion terms. The grid has been refined 
until numerical convergence was observed. 

The thick black symbols mark the benchmark solution at times t = Is, t = 3.16s and 
t= 10 s. The other curves are explained in the legend of each plot. 

As we can see in Figure [2(a) | and Figure [2(b) | the Pjv methods as well as the Djv ap- 
proaches lead to solutions that converge for increasing order iV to the benchmark solution. 
But comparing the order of the method that is necessary to reach a specific accuracy shows 
that the Dn approach leads to similar results with less computational effort. Additionally, 
we see that for small times (t — 1) both methods perform similarly well. But for large times 
(t = 10) the D\ solution agrees already very well with the benchmark solution while the Pi 
solution is much further away, especially in the region of the central peak. The solutions of 
order 5 in the Dn and Pjv approach (not shown) are almost identical and differ only in a few 
regions from the benchmark solution. 
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(c) Comparison: D\ vs. Pi approximations. 




(b) Djv approximations. 




(d) Comparison: D3 vs. P3 approximations. 



Figure 2: Energy distribution at time t = 1 s, t = 3.16 s and t = 10 s for P/v and Dn approxima- 
tions of different order. 
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Figure 3: Gray regions and the center area arc highly absorbing while white regions arc highly 
scattering. The radiation source is located in the hatched center region. 
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7.2 Lattice Problem 



This is an example with a complicated geometry, taken from [3] . We consider a checkerboard 
of highly scattering and highly absorbing regions on a lattice core. A graphical representation 
of the setting is shown in Figure [3] The white regions consist of a purely scattering material 
with k = Ocm -1 and a = 1cm -1 . The eleven gray regions and the central region are 
purely absorbing with k = 10 cm" 1 and a = Ocm -1 . For the propagation speed we assume 
c — 1 cm/s. At time zero, a source of strength one is turned on in the hatched central region. 
The computational domain is surrounded on all sides by vacuum boundaries. 

The numerical results presented here have been obtained using a finite element discretiza- 
tion with streamline diffusion. We used between 25000 and 400000 bilinear elements. More 
details on the method can be found in [21] . 

In Figure [4] we present the energy distribution ( J s2 /(fi)dfi) of the radiative field 3.2 
seconds after the radiation source in the center is turned on. The scale is logarithmic (log 10 ). 
We compare Pjv methods with Dn methods of different order. The main differences in the 
solutions can be found in the beams leaking between the corners of the absorbing regions, 
the shadows behind the absorbing regions and the front of photons escaping from the source 
region. 

As we can see in the resulting figures, for increasing order, both approaches converge to 
the same solution which for the Pr and D? models is almost the same as the one obtained 
by Monte Carlo simulations in [3]. But the Dn model gives much better results for lower 
order approximations than the Pn model. In particular, the front of the escaping photons is 
tracked much better and the shadows behind the absorbing regions are more visible in lower 
order Dn approaches. 

The fact that the front of photons is not captured that well by the Pjv method is related 
to the hyperbolic structure of the equations. Especially in the Pi model the information can 
be distributed only with one characteristic speed of l/v3 cm/s. But that is far too slow, 
compared to the real speed of the photons (lcm/ s). The higher the order N of the Pn 
approximations the more the characteristic speed of the equations approaches the desired one 
and therefore the front can be tracked much better (see [25]). Due to the additional diffusive 
term introduced by the deviation approximation into the Dm approximation, this effect is not 
present there. 
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A Properties of Spherical Harmonics 

Spherical harmonics are used as a set of basis functions for representing functions mapping 
from the unit sphere into the complex numbers. A good overview on the basics and properties 
of spherical harmonics can be found in [2]. Here only the most important properties related 
to moment methods will be recalled. 

Since Spherical harmonics act on the unit sphere, a parameterization is needed 



S 2 = {!lel 3 : fi= ^T~^sanp , ipe [0,2tt], fie [-1,1]} . (A.l) 

Using these coordinates and the associated Legendre polynomials gives rise to 
Definition A.l. For all n £ No and m £ {— n, . . . , n} the function 

Y™ : S 2 -> C 



{tp,p) i-> (-1) «/— j— - — -P„ (n e 

y 4-7T (n + my. 

is called a spherical harmonic function of order n and degree m, where the associated Legendre 
polynomials of order n and degree m are defined from the Legendre polynomials P n as 

P^{^) = {l-^—P n ^) n £ No, m — 0, . . . ,n . (A.3) 

If it is clear from the context, the dependence on Q will be neglected and we write Y™ = 
Y™(Sl). For indices m ^ {— n, . . . , n} the spherical harmonics are identically zero. For the set 
of all spherical harmonics we introduce 

Definition A.2. The vector of all spherical harmonics is given by 

Y(n) = (y (n),yr 1 (fi),yi°(fi),n 1 («),...) T c/ 2 (L 2 (5 2 ,c)) . (A.4) 

The spherical harmonics up to order N can be represented by 

Y N (n) =(Y °(tl),Yr 1 (n),Y 1 °(n),...,Y£- 1 (CI), Y$((l)y . (A.5) 

Again, if it is clear what we are referring to we neglect the dependence on Q and simply 
write Y and Yjv. 

We use the following properties of spherical harmonics to derive moment models. There 
is a relation between spherical harmonics and their complex conjugated counterpart 



Y„ m (fi) = (-l) m Y-" l (fi) , (A.6) 

an addition theorem which leads to a relation between spherical harmonics and Legendre 
polynomials 

p n (n ■ n') = V Y*(n) WJrj , (a.7) 

Zn + 1 z — ' 
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and the recursion relations 



e - Iv sin6>y r 



n 



hi(m,n)Y^ +1 1 - h(m,n)Y^ n _ 1 1 



sin# Y™ = — h2(m, n)Y^ + l2(m,n)Y^X 
cos6> Y™ = hz(m,n)Y™ +1 + h(m,n)Y^ l _ 1 



(A.8) 



with the coefficients 
hi(m, n) - 

h.2{m,n) - 

h 3 (m,n) -- 



(n — m + 1) (n — m + 2) 
(2n + l)(2n + 3) 

(n + m + l)(n + m + 2) 
(2n + l)(2n + 3) 

(n — m + l)(n + m + 1) 



li(m, n) 
Z 2 (m, n) 
h{m, n) 



(2n + l)(2?i + 3) 
It is easy to see that for the coefficients in (|A,9[l it holds 



(n + m) (n + m — 1) 
(2n- l)(2n + l) 

(n — m) (n — m — 1) 
(2n- l)(2n + l) 

(n — m)(n + m) 
(2n- l)(2n + l) ' 



(A.9) 



hi(m,n) — h2{—m,n) 
h 3 {m,n) = hz{—m,n) 



h(m,n) — l2{—m,n) 
l 3 (m,n) = l 3 (—m,n) 



(A.10) 



Using the given recursion relations leads for any vector Q on the unit sphere to 



\ (/u(™,«)^+i +/ l2 (m,n)y-+ 1 -h{ra,n)Y^i -h{m,n)Y^+ 1 ) . (A.ll) 



To express the relations (|A.11[) as matrix-vector multiplications we introduce some new 



matrices. For j £ {1,2,3} the matrices L x . £ 
defined as 



and EC 



and 



(r,s) 



(r,s) 



(r,s) 



(r,s) 



(r,s) 



(r,s) 



—li(r — k — l,k) for r = s + 2 
Z 2 ( J" — k — 1 , k) for r = s 
otherwise 

—h {r — k — 1, k) for r = s + 2 
-Z2(r — k — 1, fc) for r = s 
otherwise 

Za(r — fe — 1, k) for r = s + 1 
otherwise 



hi(r — k, k — 1) for r — s 
— fi2(r — k, k — 1) for r = s — 2 
otherwise 

hi(r — k, k — 1) for r = s 

Ii2(r — k, k — 1) for r = s — 2 , 

^ otherwise 

h 3 (r — k, k — 1) for r = s — 1 

otherwise 



(A. 12a) 

(A.12b) 
(A.12c) 

(A. 13a) 

(A.13b) 
(A.13c) 



We used the definitions of U{-, ■) and /].;(•, •) from ()A.9|I . The resulting matrices have only two 
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(or for L* and [/* only one) diagonals that do not vanish. Their structure is 



U. 



(hi 
hi 

\ 



-h 2 




ut = 



IK 



/0 h 3 

h z 



/hi h 2 

hi h 2 

V 



(A. 14) 



V 



L X1 — 



(h 












f° 




\ 


«2 






-Z 2 






h 







-h 






-h 




T k - 

U x 3 — 





h 




-il 






-Ii 













V 


'■■) 






'■■) 




I 




■J 



(A.15) 



In this representation we neglected the dependence of the coefficients li , hi , . . . on the order 
of the spherical harmonics. These relations can be found in (fA~T2)) and (|A~13l) . 

The matrices can be combined to larger matrices L^. and t) x . in the following way 



/ 



\ 







V 



/o f/i 



and tXf 



Additionally we introduce the factor 



7i 



\ 



X 3 

o / 



(A.16) 



for j — 1 
for j = 2 
1 for j = 3 



and the unit vector 



9ejl ! „ ) = (0,...,0,l,0,...,0) r 



(A. 17) 



(A.18) 



which is one only in the component that is related to the spherical harmonic Y™ in the vector 



Yat such that (e 



(m,n) 



write 



= Y™. Then, for n € {0, . . . , N} and m G {— n, . . . , n}, we can 
r= 7 ,((ef m , n) ) T ^.Y^+(e^ ) ) T ^ +1 Y JV+1 ) . (A.19) 



Lemma A. 3. For the matrices L T . and U. r . hold the relations 



(A.20) 



B Treatment of Scattering and Absorption opera- 
tors in moment methods 



The radiative transfer equation contains a scattering component 

(si)(t,x,n) = — / $(x,n-n')i(t,x.n')drt (b.i) 

47T Jg2 

with a normalized scattering kernel <!>. This scattering kernel can be rather complicated. Sev- 
eral theories have been developed to approximate realistic kernels. The first who established 
an approach was Mie [24| . However, due to the complexity of his theory, several simplified 
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approximations have been developed, e.g. the Henyey-Greenstein (HG) kernel [28] or the SAM 
(simplified approximate Mie) approach [16] . 

To deal with the scattering kernel in moment methods, it is very common to rewrite it 
into a series expansion based on Legendre polynomials 



2 

3=0 

This leads to 

oo I 



*(x,ri=j£y±±<Ti(x)PM- (B.2) 

oo I 

5/(0) = |£ Y, °tfY?(Sl) (B.3) 



and 

OO I OO I 



m r a SI(Q) = ?EE ("WO)) = ? E E = f • (B-4) 



i=0 fc=-( i=0 *=-! 



For the total absorption and scattering operator we get 

(O)(o) = E ( K + a - f °0 ^^(n) = : E E (b.5) 



oo i oo Z 

(7 

2' 

i=0 fc=— J i=0 fc=— 1 

It is obvious that the inverse operator is 

oo l 

(/c- 1 /)^) = e E t 5 i<y» fc (") • (b.6) 

and thus 

mI(/C- 1 /)(n) = ■ 1 - II. (B.7) 
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